Asaia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
path_asaia <- "../../../../output/2_Oligotyping/2A/Asaia/2A_oligotyping_Asaia_sequences-c5-s1-a0.0-A0-M10"

# load the matrix count table
matrix_count <- read.table(paste0(path_asaia, "/MATRIX-COUNT.txt"), header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP14 NP2 NP20 NP27 NP29 NP30 NP34 NP35 NP36 NP37 NP38 NP39 NP41 NP42 NP43 NP44 NP5 NP8 S126 S147 S154 S18 S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 S42 S44 S45 S46 S47 S48 S49 S50 S51 S52
TGCGA 0 1 1 0 8 0 0 0 1 1 0 2 0 0 32 0 70 0 73 70 11 8 221 153738 118463 81849 97920 46367 64795 135582 0 0 0 0 0 0 2 0 0 0 0 1 1 0 6 0 0 0 1 1 0 0 0 7 10 5 15 18 0 8 0 0 0 0 0 32 0 0
CATAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 31 0 0 1 58 56 8 7 2 2744 129509 116951 141940 121248 146560 1893 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TACGA 57 315 647 156 2206 18 133 9 567 377 1 2165 12 88 6 390 1 5 43 13 8 11 7 196884 20230 16952 3701 15217 22358 17187 1 16 1 1 2 17 1277 67 195 21 33 468 108 109 866 105 1 1006 1355 1772 2848 1817 5507 6828 9095 3991 18210 20197 1 24 46 24 84 85 105 182 310 298
TGCGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 4 0 0 2 718 627 479 857 452 679 1193 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CATGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 35 779 731 1039 900 1113 21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TACGG 0 2 0 0 5 0 0 0 1 1 0 7 0 0 0 3 0 0 1 7 0 0 0 844 175 112 72 153 308 157 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 2 0 2 0 2 7 2 4 26 20 0 0 0 0 0 0 0 0 0 0
CATAA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 25 72 52 59 70 1 0 0 0 0 0 0 7 2 3 0 1 10 2 4 3 4 0 23 55 93 68 27 193 55 87 31 92 359 0 0 1 0 1 1 3 0 0 5
CGTAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 1 99 138 169 134 242 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CGCGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 74 78 111 89 29 53 174 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CATGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 35 33 135 45 44 61 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 2 0 0 0 0 0 0 0 0 0 0
TGCAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 17 31 116 64 36 55 17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CGTGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 4 57 33 67 74 82 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TATAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 11 66 59 56 60 63 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TGTGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 67 62 27 41 22 44 45 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TGCAA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 81 26 26 38 20 23 63 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TATGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 21 45 27 33 51 35 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CACAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 16 41 38 37 38 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TATGA 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 43 13 12 20 18 16 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 4 0 0 0 6 4 4 6 4 0 0 0 0 0 0 0 0 0 0
CACGA 0 0 0 0 3 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 62 19 7 4 5 14 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 1 0 0 0 0 0 2 4 0 15 9 0 0 0 0 0 0 0 0 0 0
TACAA 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 72 2 8 3 6 9 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 3 2 11 0 6 19 0 0 0 0 0 0 0 0 0 1
TGTGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 13 24 10 18 39 27 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TACAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 21 7 31 2 17 23 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TGCTA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 10 11 22 11 11 26 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TACTA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 22 4 5 3 0 6 3 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 2 3 1 1 4 2 6 0 10 23 0 0 0 0 0 0 1 0 0 0
CGCGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 10 6 14 14 18 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CTTAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 11 11 12 19 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CACGG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 6 5 10 7 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TGTAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 6 3 3 5 15 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
GGCGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 25 0 7 2 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CATAT 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 2 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CGTGA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 5 0 8 4 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
AATAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3 6 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TATAA 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 5 3 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CGCAG 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 2 3 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Taxonomy

# load the fasta table
fasta <- readDNAStringSet(paste0(path_asaia, "/OLIGO-REPRESENTATIVES.fasta"))

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
oligotype
TGCGA TGCGA
CATAG CATAG
TACGA TACGA
TGCGG TGCGG
CATGG CATGG
TACGG TACGG
CATAA CATAA
CGTAG CGTAG
CGCGA CGCGA
CATGA CATGA
TGCAG TGCAG
CGTGG CGTGG
TATAG TATAG
TGTGA TGTGA
TGCAA TGCAA
TATGG TATGG
CACAG CACAG
TATGA TATGA
CACGA CACGA
TACAA TACAA
TGTGG TGTGG
TACAG TACAG
TGCTA TGCTA
TACTA TACTA
CGCGG CGCGG
CTTAG CTTAG
CACGG CACGG
TGTAG TGTAG
GGCGA GGCGA
CATAT CATAT
CGTGA CGTGA
AATAG AATAG
TATAA TATAA
CGCAG CGCAG

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## load the reference table
ref_oligo_med2 <- read.table("../../../../output/2_Oligotyping/2B/2B_REF_info_asaia.tsv", sep="\t", header = TRUE)

## select only the 5 oligotypes of Asaia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP14 NP2 NP20 NP27 NP29 NP30 NP34 NP35 NP36 NP37 NP38 NP39 NP41 NP42 NP43 NP44 NP5 NP8 S126 S147 S154 S18 S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 S42 S44 S45 S46 S47 S48 S49 S50 S51 S52
oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340 0 0 0 0 0 0 0 0 0 0 0 0 0 0 31 0 0 1 58 56 8 7 2 2744 129509 116951 141940 121248 146560 1893 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
oligotype_CATGG (8) | size:4619 / node_N0492 (18) | size:6153 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 35 779 731 1039 900 1113 21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329 57 315 647 156 2206 18 133 9 567 377 1 2165 12 88 6 390 1 5 43 13 8 11 7 196884 20230 16952 3701 15217 22358 17187 1 16 1 1 2 17 1277 67 195 21 33 468 108 109 866 105 1 1006 1355 1772 2848 1817 5507 6828 9095 3991 18210 20197 1 24 46 24 84 85 105 182 310 298
oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167 0 1 1 0 8 0 0 0 1 1 0 2 0 0 32 0 70 0 73 70 11 8 221 153738 118463 81849 97920 46367 64795 135582 0 0 0 0 0 0 2 0 0 0 0 1 1 0 6 0 0 0 1 1 0 0 0 7 10 5 15 18 0 8 0 0 0 0 0 32 0 0
oligotype_TGCGG (10) | size:5016 / node_N1146 (10) | size:4924 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 4 0 0 2 718 627 479 857 452 679 1193 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
ref
oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340 oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340
oligotype_CATGG (8) | size:4619 / node_N0492 (18) | size:6153 oligotype_CATGG (8) | size:4619 / node_N0492 (18) | size:6153
oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329 oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329
oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167 oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167
oligotype_TGCGG (10) | size:5016 / node_N1146 (10) | size:4924 oligotype_TGCGG (10) | size:5016 / node_N1146 (10) | size:4924

Metadata

metadata <- read.csv("../../../../metadata/metadata.csv", sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 5 taxa and 68 samples ]
## sample_data() Sample Data:       [ 68 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 5 taxa by 1 taxonomic ranks ]
compute_read_counts(ps)
## [1] 1746802
# remove blanks
ps <- subset_samples(ps, Strain!="Blank")
ps <- check_ps(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 5 taxa and 67 samples ]
## sample_data() Sample Data:       [ 67 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 5 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_asaia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Asaia
ps.asaia <- ps.filter %>% subset_taxa(Genus=="Asaia")

# add read depth of Asaia
sample_data(ps.filter)$Read_asaia <- sample_sums(ps.asaia)
sample_data(ps.filter) %>% colnames()
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"        
## [13] "Read_asaia"
sample_data(ps.asaia) %>% colnames()
##  [1] "Sample"         "Strain"         "Field"          "Country"       
##  [5] "Organ"          "Species"        "Run"            "Control"       
##  [9] "Species_italic" "Strain_italic"  "Read_depth"     "is.neg"
# add percent of Asaia
sample_data(ps.filter)$Percent_asaia <- sample_data(ps.filter)$Read_asaia / sample_data(ps.filter)$Read_depth

# round the percent of Asaia at 2 decimals
sample_data(ps.filter)$Percent_asaia <- sample_data(ps.filter)$Percent_asaia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_asaia, Percent_asaia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Sample Strain Field Country Organ Species Run Control Species_italic Strain_italic Read_depth is.neg Read_asaia Percent_asaia
CTC1 CTC1 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 29779 FALSE 57 0.00
CTC10 CTC10 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2609 FALSE 318 0.12
CTC11 CTC11 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 13874 FALSE 649 0.05
CTC12 CTC12 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1146 FALSE 156 0.14
CTC13 CTC13 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 18035 FALSE 2224 0.12
CTC14 CTC14 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1708 FALSE 18 0.01
CTC15 CTC15 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 23180 FALSE 133 0.01
CTC2 CTC2 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 30692 FALSE 9 0.00
CTC3 CTC3 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 39920 FALSE 569 0.01
CTC4 CTC4 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2139 FALSE 380 0.18
CTC5 CTC5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 15789 FALSE 1 0.00
CTC6 CTC6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 19753 FALSE 2178 0.11
CTC9 CTC9 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus run2 True sample italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 4980 FALSE 88 0.02
NP14 NP14 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 7973 FALSE 69 0.01
NP2 NP2 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 648335 FALSE 394 0.00
NP20 NP20 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 136 FALSE 71 0.52
NP27 NP27 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 1234 FALSE 6 0.00
NP29 NP29 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 203 FALSE 184 0.91
NP30 NP30 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 228 FALSE 155 0.68
NP34 NP34 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 95 FALSE 27 0.28
NP35 NP35 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 196532 FALSE 26 0.00
NP36 NP36 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 249 FALSE 232 0.93
NP37 NP37 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 419340 FALSE 355588 0.85
NP38 NP38 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 282479 FALSE 270454 0.96
NP39 NP39 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 218684 FALSE 218013 1.00
NP41 NP41 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 247152 FALSE 246366 1.00
NP42 NP42 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 185157 FALSE 185069 1.00
NP43 NP43 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 239335 FALSE 236790 0.99
NP44 NP44 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti run3 True sample italic(“Aedes aegypti”) Field-Guadeloupe 156879 FALSE 156424 1.00
NP5 NP5 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 736159 FALSE 2 0.00
NP8 NP8 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus run3 True sample italic(“Culex quinquefasciatus”) Field-Guadeloupe 334799 FALSE 17 0.00
S100 S100 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52486 FALSE 0 0.00
S102 S102 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 3456 FALSE 0 0.00
S104 S104 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52403 FALSE 0 0.00
S105 S105 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 55577 FALSE 0 0.00
S106 S106 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 33053 FALSE 0 0.00
S107 S107 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52154 FALSE 0 0.00
S108 S108 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 55735 FALSE 0 0.00
S109 S109 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 59023 FALSE 0 0.00
S110 S110 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 57377 FALSE 0 0.00
S121 S121 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 20361 FALSE 0 0.00
S122 S122 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 9803 FALSE 0 0.00
S123 S123 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 20130 FALSE 0 0.00
S124 S124 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 18146 FALSE 0 0.00
S126 S126 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 15235 FALSE 1 0.00
S127 S127 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 24696 FALSE 0 0.00
S128 S128 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 16305 FALSE 0 0.00
S146 S146 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25012 FALSE 0 0.00
S147 S147 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25171 FALSE 1 0.00
S148 S148 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 14164 FALSE 0 0.00
S150 S150 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 15081 FALSE 0 0.00
S151 S151 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 22944 FALSE 0 0.00
S152 S152 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 15082 FALSE 0 0.00
S153 S153 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 17040 FALSE 0 0.00
S154 S154 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 9626 FALSE 2 0.00
S160 S160 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 72508 FALSE 0 0.00
S162 S162 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 25180 FALSE 0 0.00
S163 S163 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 12333 FALSE 0 0.00
S164 S164 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 22368 FALSE 0 0.00
S165 S165 Laboratory - Lavar Lab France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Laboratory-Lavar 17731 FALSE 0 0.00
S166 S166 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 13979 FALSE 0 0.00
S167 S167 Field - Bosc Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Bosc 14048 FALSE 0 0.00
S169 S169 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 11553 FALSE 0 0.00
S170 S170 Field - Camping Europe Field France Ovary Culex pipiens run2 True sample italic(“Culex pipiens”) Field-Camping~Europe 8852 FALSE 0 0.00
S18 S18 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 4290 FALSE 17 0.00
S19 S19 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 44527 FALSE 1288 0.03
S20 S20 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 42864 FALSE 69 0.00
S21 S21 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33798 FALSE 198 0.01
S22 S22 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 19044 FALSE 21 0.00
S23 S23 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 38172 FALSE 34 0.00
S24 S24 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 42355 FALSE 479 0.01
S25 S25 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 47688 FALSE 111 0.00
S26 S26 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 5394 FALSE 113 0.02
S27 S27 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 24558 FALSE 879 0.04
S28 S28 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 4503 FALSE 109 0.02
S30 S30 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 25353 FALSE 1032 0.04
S31 S31 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 20417 FALSE 1416 0.07
S32 S32 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 12441 FALSE 1873 0.15
S33 S33 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33867 FALSE 2919 0.09
S34 S34 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 9367 FALSE 1845 0.20
S35 S35 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 11663 FALSE 5709 0.49
S36 S36 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 33020 FALSE 6912 0.21
S37 S37 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 18340 FALSE 9220 0.50
S38 S38 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 54790 FALSE 4036 0.07
S39 S39 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 36273 FALSE 18381 0.51
S40 S40 Laboratory - Lavar Lab France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Laboratory-Lavar 44448 FALSE 20654 0.46
S42 S42 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 4107 FALSE 1 0.00
S43 S43 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 9279 FALSE 0 0.00
S44 S44 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 8026 FALSE 32 0.00
S45 S45 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 18150 FALSE 47 0.00
S47 S47 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 1951 FALSE 85 0.04
S48 S48 Field - Camping Europe Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 56738 FALSE 86 0.00
S49 S49 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 33498 FALSE 109 0.00
S50 S50 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 28481 FALSE 214 0.01
S51 S51 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 61788 FALSE 310 0.01
S52 S52 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 21553 FALSE 304 0.01
S55 S55 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 50447 FALSE 0 0.00
S56 S56 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 42609 FALSE 0 0.00
S57 S57 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 49157 FALSE 0 0.00
S58 S58 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 30357 FALSE 0 0.00
S59 S59 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 32798 FALSE 0 0.00
S60 S60 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 44485 FALSE 0 0.00
S61 S61 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 49545 FALSE 0 0.00
S63 S63 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 53444 FALSE 0 0.00
S64 S64 Field - Bosc Field France Whole Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 47628 FALSE 0 0.00
S79 S79 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 59755 FALSE 0 0.00
S80 S80 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 52788 FALSE 0 0.00
S83 S83 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 42272 FALSE 0 0.00
S84 S84 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 56676 FALSE 0 0.00
S85 S85 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 41690 FALSE 0 0.00
S86 S86 Field - Camping Europe Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Camping~Europe 61984 FALSE 0 0.00
S87 S87 Field - Bosc Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 65958 FALSE 0 0.00
S88 S88 Field - Bosc Field France Ovary Culex pipiens run1 True sample italic(“Culex pipiens”) Field-Bosc 53102 FALSE 0 0.00 
# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Percent by oligotype

# Guadeloupe 

## Oligotype TGCGA

oligo_TGCGA <- ps %>% 
  subset_samples(Organ=="Whole" & Strain=="Field - Guadeloupe") %>%
  subset_taxa(ref=="oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167")

oligo_TGCGA <- prune_taxa(taxa_sums(oligo_TGCGA) >= 1, oligo_TGCGA)
oligo_TGCGA <- prune_samples(sample_sums(oligo_TGCGA) >= 1, oligo_TGCGA)
  
oligo_TGCGA %>% taxa_sums() -> sum_oligo_1
oligo_TGCGA@sam_data$Read_depth %>% sum() -> sum_total_1

sum_oligo_1 / sum_total_1 *100 # 35.91%
## oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167 
##                                                           35.91867
sum_oligo_1 / 6452623 *100 # 10.83%
## oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167 
##                                                           10.83431
## Oligotype CATAG

oligo_CATAG <- ps %>% 
  subset_samples(Organ=="Whole" & Strain=="Field - Guadeloupe") %>%
  subset_taxa(ref=="oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340")

oligo_CATAG <- prune_taxa(taxa_sums(oligo_CATAG) >= 1, oligo_CATAG)
oligo_CATAG <- prune_samples(sample_sums(oligo_CATAG) >= 1, oligo_CATAG)

oligo_CATAG %>% taxa_sums() -> sum_oligo_2
oligo_CATAG@sam_data$Read_depth %>% sum() -> sum_total_2

sum_oligo_2 / sum_total_2 *100 # 33.94%
## oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340 
##                                                            33.9386
sum_oligo_2 / 6452623 *100 # 10.24%
## oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340 
##                                                           10.24354
## Oligotype TACGA

oligo_TACGA <- ps %>% 
  subset_samples(Organ=="Whole" & Strain=="Field - Guadeloupe") %>%
  subset_taxa(ref=="oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329")

oligo_TACGA <- prune_taxa(taxa_sums(oligo_TACGA) >= 1, oligo_TACGA)
oligo_TACGA <- prune_samples(sample_sums(oligo_TACGA) >= 1, oligo_TACGA)
  
oligo_TACGA %>% taxa_sums() -> sum_oligo_3 
oligo_TACGA@sam_data$Read_depth %>% sum() -> sum_total_3 

sum_oligo_3 / sum_total_3 *100 # 15.02%
## oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329 
##                                                           15.02469
sum_oligo_3 / 6452623 *100 # 4.53%
## oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329 
##                                                           4.534838
# France

## Oligotype TGCGA

oligo_TGCGA <- ps %>% 
  subset_samples(Organ=="Whole" & Strain!="Field - Guadeloupe") %>%
  subset_taxa(ref=="oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167")

oligo_TGCGA <- prune_taxa(taxa_sums(oligo_TGCGA) >= 1, oligo_TGCGA)
oligo_TGCGA <- prune_samples(sample_sums(oligo_TGCGA) >= 1, oligo_TGCGA)
  
oligo_TGCGA %>% taxa_sums() -> sum_oligo_1
oligo_TGCGA@sam_data$Read_depth %>% sum() -> sum_total_1

sum_oligo_1 / sum_total_1 *100 # 0.024%
## oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167 
##                                                         0.02364695
sum_oligo_1 / 6452623 *100 # 0.0019%
## oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167 
##                                                        0.001875206
## Oligotype CATAG

oligo_CATAG <- ps %>% 
  subset_samples(Organ=="Whole" & Strain!="Field - Guadeloupe") %>%
  subset_taxa(ref=="oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340")

oligo_CATAG <- prune_taxa(taxa_sums(oligo_CATAG) >= 1, oligo_CATAG)
oligo_CATAG <- prune_samples(sample_sums(oligo_CATAG) >= 1, oligo_CATAG)

oligo_CATAG %>% taxa_sums() -> sum_oligo_2
oligo_CATAG@sam_data$Read_depth %>% sum() -> sum_total_2

sum_oligo_2 / sum_total_2 *100 # 0.003%
## oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340 
##                                                        0.003028468
sum_oligo_2 / 6452623 *100 # 1.5e-05 %
## oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340 
##                                                       1.549757e-05
## Oligotype TACGA

oligo_TACGA <- ps %>% 
  subset_samples(Organ=="Whole" & Strain!="Field - Guadeloupe") %>%
  subset_taxa(ref=="oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329")

oligo_TACGA <- prune_taxa(taxa_sums(oligo_TACGA) >= 1, oligo_TACGA)
oligo_TACGA <- prune_samples(sample_sums(oligo_TACGA) >= 1, oligo_TACGA)
  
oligo_TACGA %>% taxa_sums() -> sum_oligo_3 
oligo_TACGA@sam_data$Read_depth %>% sum() -> sum_total_3 

sum_oligo_3 / sum_total_3 *100 # 8.01%
## oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329 
##                                                           8.015364
sum_oligo_3 / 6452623 *100 # 1.3%
## oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329 
##                                                            1.29817

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Strain, Organ, Read_depth, Read_asaia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Strain", "Organ"), vars=c("Read_depth", "Read_asaia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
# afficher plot
p_count1
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0), nrow=3, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 5 taxa and 57 samples ]
## sample_data() Sample Data:       [ 57 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 5 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
data_for_plot2 <- taxo_data(ps.filter.whole, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()
## "oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340.",
##  "oligotype_CATGG (8) | size:4619 / node_N0492 (18) | size:6153.",
##  "oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329.",
##  "oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167.",
##  "oligotype_TGCGG (10) | size:5016 / node_N1146 (10) | size:4924.",
data_for_plot2$Name <- data_for_plot2$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot2$Name <- as.factor(data_for_plot2$Name)

new_names <- c("oligotype_TGCGA (33) | size:699320 / N1147 (33) | size:651167.",
               "oligotype_CATAG (16) | size:661010 / N1156 (40) | size:622340.",
               "oligotype_TACGA (68) | size:376837 / N0939 (68) | size:364329.",
               "oligotype_TGCGG (10) | size:5016 / N1146 (10) | size:4924.",
               "oligotype_CATGG (8) | size:4619 / N0492 (18) | size:6153.",
               "Other.")


data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col_add <- brewer.pal(8, "Accent")


col <- c("oligotype_TGCGA (33) | size:699320 / N1147 (33) | size:651167."="#CDDE8E",
         "oligotype_CATAG (16) | size:661010 / N1156 (40) | size:622340."="#DDFFC4",
         "oligotype_TACGA (68) | size:376837 / N0939 (68) | size:364329."="#A6DE45",
         "oligotype_TGCGG (10) | size:5016 / N1146 (10) | size:4924."=col_add[1],
         "oligotype_CATGG (8) | size:4619 / N0492 (18) | size:6153."=col_add[2],
         "Other."="#A0A0A0")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Asaia"), "- (Slab TC)")))

#data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("")+
  guide_italics+
  theme(legend.title = element_text(size = 20),
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"),
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3 taxa and 8 samples ]
## sample_data() Sample Data:       [ 8 samples by 14 sample variables ]
## tax_table()   Taxonomy Table:    [ 3 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot3 <- taxo_data_fast(ps.filter.ovary, method = "abundance")
data_for_plot3 <- taxo_data(ps.filter.ovary, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()
## "oligotype_CATAG (16) | size:661010 / node_N1156 (40) | size:622340.",
##  "oligotype_TACGA (68) | size:376837 / node_N0939 (68) | size:364329.",
##  "oligotype_TGCGA (33) | size:699320 / node_N1147 (33) | size:651167.",
data_for_plot3$Name <- data_for_plot3$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot3$Name <- as.factor(data_for_plot3$Name)

new_names <- c("oligotype_TGCGA (33) | size:699320 / N1147 (33) | size:651167.",
               "oligotype_CATAG (16) | size:661010 / N1156 (40) | size:622340.",
               "oligotype_TACGA (68) | size:376837 / N0939 (68) | size:364329.",
               "Other."
)

data_for_plot3$Name <- factor(data_for_plot3$Name, levels = new_names)

levels(data_for_plot3$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

#levels(data_for_plot3$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Asaia"), "- (Slab TC)")))
levels(data_for_plot3$Strain) <- c("Bosc", "Guadeloupe", "Lavar~(lab)")

#data_for_plot3 <- data_for_plot3 %>% na.omit()

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

tiff("../../../../output/2_Oligotyping/2D/2Db_OLIGO_counts_asaia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_counts_asaia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_counts_asaia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_taxonomic_asaia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2

Make main plot

img <- magick::image_read(paste0(path_asaia, "/HTML-OUTPUT/entropy.png"))
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

tiff("../../../../output/2_Oligotyping/2D/2Db_OLIGO_main_asaia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/2_Oligotyping/2D/2Db_OLIGO_main_asaia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2

Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.1    openxlsx_4.2.3      Biostrings_2.54.0  
##  [4] XVector_0.26.0      IRanges_2.20.2      S4Vectors_0.24.4   
##  [7] BiocGenerics_0.32.0 cowplot_1.1.0       ggpubr_0.4.0       
## [10] RColorBrewer_1.1-2  forcats_0.5.0       stringr_1.4.0      
## [13] dplyr_1.0.2         purrr_0.3.4         readr_1.4.0        
## [16] tidyr_1.1.2         tibble_3.0.4        tidyverse_1.3.0    
## [19] ggplot2_3.3.2       phyloseq_1.30.0    
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-0  ggsignif_0.6.0    ellipsis_0.3.1    rio_0.5.16       
##  [5] fs_1.5.0          rstudioapi_0.11   farver_2.0.3      fansi_0.4.1      
##  [9] lubridate_1.7.9   xml2_1.3.2        codetools_0.2-16  splines_3.6.3    
## [13] knitr_1.30        ade4_1.7-15       jsonlite_1.7.1    broom_0.7.2      
## [17] cluster_2.1.0     dbplyr_1.4.4      compiler_3.6.3    httr_1.4.2       
## [21] backports_1.1.10  assertthat_0.2.1  Matrix_1.2-18     cli_2.1.0        
## [25] htmltools_0.5.1.1 tools_3.6.3       igraph_1.2.6      gtable_0.3.0     
## [29] glue_1.4.2        reshape2_1.4.4    Rcpp_1.0.5        carData_3.0-4    
## [33] Biobase_2.46.0    cellranger_1.1.0  jquerylib_0.1.3   vctrs_0.3.4      
## [37] multtest_2.42.0   ape_5.4-1         nlme_3.1-149      iterators_1.0.13 
## [41] xfun_0.22         ps_1.4.0          rvest_0.3.6       lifecycle_0.2.0  
## [45] rstatix_0.6.0     zlibbioc_1.32.0   MASS_7.3-53       scales_1.1.1     
## [49] hms_0.5.3         biomformat_1.14.0 rhdf5_2.30.1      yaml_2.2.1       
## [53] curl_4.3          sass_0.3.1        stringi_1.5.3     highr_0.8        
## [57] foreach_1.5.1     permute_0.9-5     zip_2.1.1         rlang_0.4.10     
## [61] pkgconfig_2.0.3   evaluate_0.14     lattice_0.20-41   Rhdf5lib_1.8.0   
## [65] labeling_0.4.2    tidyselect_1.1.0  plyr_1.8.6        magrittr_1.5     
## [69] bookdown_0.22     R6_2.4.1          magick_2.5.2      generics_0.0.2   
## [73] DBI_1.1.0         pillar_1.4.6      haven_2.3.1       foreign_0.8-75   
## [77] withr_2.3.0       mgcv_1.8-33       survival_3.2-7    abind_1.4-5      
## [81] modelr_0.1.8      crayon_1.3.4      car_3.0-10        rmarkdown_2.7    
## [85] grid_3.6.3        readxl_1.3.1      data.table_1.13.2 blob_1.2.1       
## [89] vegan_2.5-6       rmdformats_1.0.2  reprex_0.3.0      digest_0.6.26    
## [93] webshot_0.5.2     munsell_0.5.0     viridisLite_0.3.0 bslib_0.2.4